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Abstract 

One of the most interesting questions regarding a possible first order 
cosmological quark-hadron phase transition concerns the final fate of 
the baryon number contained within the disconnected quark regions 
at the end of the transition. We here present a detailed investigation 
of the hydro dynamical evolution of an evaporating quark drop, using a 
multi-component fluid description to follow the mechanisms of baryon 
number segregation. With this approach, we are able to take account 
of the simultaneous effects of baryon number flux suppression at the 
phase interface, entropy extraction by means of particles having long 
mean- free-paths, and baryon number diffusion. A range of computa- 
tions has been performed to investigate the permitted parameter-space 
and this has shown that significant baryon number concentrations, 
perhaps even up to densities above that of nuclear matter, represent 
an inevitable outcome within this scenario. 
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I. Introduction 



If the cosmological quark-hadron phase transition was of first order with non- 
neghgible supercoohng (as we will be assuming in this paper), then it almost cer- 
tainly produced inhomogeneities of some kind in the spatial distribution of baryon 
number. These inhomogeneities could have had a relevant influence on the subse- 
quent nucleosynthesis if they were produced on large enough scales, had sufficiently 
large amplitude and contained a significant fraction of the baryon number present 
in the Universe 

Many authors in the past years have calculated the possible consequences 
of baryon number inhomogeneities for the subsequent evolution of the Universe, 
focussing particularly on the scenarios for inhomogeneous nucleosynthesis ^. 
However, in all of these studies it was necessary to introduce suitably chosen pa- 
rameters in order to compensate for lack of knowledge about the spatial distribution 
and amplitude of the baryon number peaks which might have been left behind at 
the end of the quark-hadron transition. There are various origins for this lack of 
knowledge which we will now discuss briefiy. A first and major unresolved issue 
concerns the typical distance between nucleation sites at the beginning of the tran- 
sition. This length scale is not only relevant at the time when hadron bubbles 
are nucleated within the slightly supercooled quark medium, but also represents a 
characteristic length scale for subsequent stages of the transition. It is the scale 
at which neighbouring hadron bubbles collide and percolate and it also determines 
the typical distance between the centres of the disconnected quark regions produced 
by bubble coalescence. As a consequence of this, it also represents the maximum 
length-scale for the production of peaks in the baryon number density at the end of 
the transition. 

Independently of the scale over which baryon number fiuctuations are pro- 
duced, the study of baryon number segregation is complicated by the presence of 
various mechanisms which could contribute to the scale of the fiuctuations. When 
there is chemical equilibrium between the two phases, baryon number density is 
already "naturally" higher in the quark phase than in the hadron phase. This dif- 
ference can be further increased as a result of various other processes which probably 
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act together. One of these is suppression of baryon number flow across the phase 
interface [|10], ^ Simple statistical considerations suggest that baryon number 
cannot be carried entirely together with the hydrodynamical flow when this moves 
from the quark to the hadron phase. A phenomenological explanation for this can 
be found in the fact that it is generally more difficult to find in a volume of 1 fm^ 
and in a time of 10~^^ sec the right triplet of up and down quarks necessary to 
form a color singlet nucleon (either a baryon or an antibaryon) than it is to find the 
doublet of quarks necessary to form the lightest hadrons (e.g. pions). 

Another mechanism, which could affect the variations in baryon number den- 
sity and which is effective during the final stages of the transition, is long-range 
energy and momentum transfer by means of particles having long mean free path 
(which can be viewed as a form of radiative transfer). This process, which has re- 
cently been investigated by Rezzolla and Miller |Tl|, involves the relativistic fiuids 
of strongly interacting particles (i.e. the hadron fiuid and the quark fiuid which we 
will refer to as "standard fluids") and the relativistic fluid of particles having only 
electromagnetic and weak interactions (which we refer to as the "radiation fluid") 
with entropy being extracted from the quark phase by the long mean-free-path par- 
ticles of the radiation fluid. A necessary condition for this process to be signiflcant 
is that the dimension Rg of the hadron phase (in the case of bubble growth) or of the 
quark phase (in the case of drop evaporation) needs to be roughly comparable with 
the smallest mean free path A of the radiation fluid particlesi. When this condition 
is not satisfled, it is a good approximation to consider the two fluids as either being 
completely coupled (for Rg ^ A) or completely decoupled (for Rg <^ A) and the 
energy and momentum interchange between the standard fluid and the radiation 
fluid is then either maximally efficient or nonexistent. 

Besides the baryon number segregating mechanisms mentioned above, there is 
also another competing process which needs to be taken into account. The formation 

^Note that in general there are two of these length scales, referring to the electromag- 
netic interaction (relevant for photons, e^, etc., for which A ~ 10'* fm) and to the 
weak interaction (relevant for neutrinos and antineutrinos, for which A ~ 10^^ fm). 



2 



of localized regions of high baryon number density is counteracted by baryon number 
diffusion which occurs whenever a deviation from homogeneity is produced and is 
effective both during the quark-hadron transition and after it. 

While, from the above discussion, the formation of peaks in baryon number 
density appears to be a rather inevitable consequence of a first order quark-hadron 
transition, a number of quantitative aspects of this picture still remain to be fully 
investigated. For this reason, and before the microphysics of baryon number flow 
across the phase interface or of baryon number diffusion is known in more detail, 
we here present results from systematic numerical computations of baryon number 
concentration during the evaporation of an isolated quark drop. These simulations 
are an extension of earlier calculations of the final stages of the quark-hadron tran- 
sition, in which the hydrodynamics of a spherical isolated quark drop was coupled 
to the processes of relativistic radiative transfer between the standard fluids and 
the radiation fluid. In the present simulations, we have also studied the behaviour 
of the additional distinct "fluid" of baryon number carriers and have followed the 
progress of the segregation as it is promoted by the suppression mechanism at the 
interface and by the entropy extraction but moderated by baryon diffusion. 

Results from computations performed with different values of the relevant 
parameters are presented in Section V. Before this. Section II contains a brief review 
of the mechanisms responsible for baryon number segregation and of their role during 
the various stages of the transition. Section III presents the system of relativistic 
equations which we use for calculating the hydrodjTiamics of the quark drop and 
the radiative transfer processes. (The formal derivation of these equations is not 
given here but can be found in references [^, 13, 14, 11 1.) In Section IV we derive a 
relativistic equation for the diffusion of baryon number and write it in a convenient 
finite difference form. This Section also introduces a phenomenological analytical 
expression by means of which the suppression mechanisms at the interface can be 
represented. Section VI contains the conclusions and also a discussion of our results 
in relation to the scenarios for inhomogeneous cosmological nucleosynthesis and for 
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the production of primordial magnetic fields. We adopt units for which c — h — 
kg — 1 and use commas or the usual "9" notation to denote partial derivatives. 

II. Evolution of the baryon number concentra- 
tion 

As mentioned in the Introduction, production of inhomogeneities in baryon 
number density can be affected by a number of different factors. In this Section 
we outline a schematic picture of the various stages of the transition during which 
baryon number segregation can either be favoured or impeded. First, we note that 
even in the absence of specific mechanisms for this segregation, there is a "natural" 
tendency for the baryon number density to be higher in the quark phase than in 
the hadron phase since baryon number is carried by almost massless quarks in the 
quark phase, whereas in the hadron phase it is carried by heavy nucleons whose 
number density is suppressed by an exponential factor. An estimate for the contrast 
in baryon number density resulting from this can easily be calculated for stages at 
which global chemical equilibrium is near to holding. The chemical potentials in 
the two phases can then be set equal {fil = = /X5) signifying that, if the phase 
interface were not in motion, equal fluxes of baryon number would cross it in both 
directions. 

Because of the complexity of QCD and of the spectrum of hadronic species 
present in the low temperature (hadron) phase, it is necessary to make approxi- 
mations if we want to obtain simple analytical expressions for the baryon number 
density in the two phases. Following the approach of our earlier papers, we use the 
phenomenological bag model for the quark phase (treating the quark-gluon plasma 
as consisting of relativistic non-interacting particles within a false vacuum) and the 
energy density and pressure may then be written as 




(1) 
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where g* is the effective number of degrees of freedom in the quark phase (including 
also the relevant number of degrees of freedom for the particles of the radiation fluid 
when the radiation and standard fluids are fully coupled [[l^])- The baryon number 
density is then given by 



NcNf 3 



h 27 VT 



T3, (2) 



where Nc and Nf are, respectively, the number of colours (three) and the number 
of quark flavours (which we will take here to be two). On a similar level of approxi- 
mation, we can consider the pion contribution to the grand partition function of the 
hadron phase as being the dominant one coming from the whole spectrum of hadron 
species, neglect finite volume effects and the repulsive interaction between hadrons 
and describe the hadron phase as being essentially a gas of relativistic massless 
particles with the equation of state 

eh = 9*h , Ph = , (3) 

where is the effective number of degrees of freedom in the hadron phase (equivalent 
to g* in the quark phase). It is then possible to obtain a simple and compact 
expression for the baryon number density in the hadron phase as 



„J . (^^f (^) e-"'- . {^f (^) . (4) 

with m being the nucleon mass. 

With all of these assumptions, the contrast in baryon number density between 
the two phases is given by 

^^!|^^/^A'^%Wt_10 for T = T, = 150 MeV, (5) 
9 \2mJ ' ^ ^ 
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where T^. is the critical temperature for the transition. Similar estimates of the 
baryon number contrast have been found also from more detailed analyses of the 
baryon number density in the two phases . 



It is important to stress that when the hypothesis of global chemical equilib- 
rium is valid, equation holds between any two generic points in the two phases. 
It is reasonable to assume that global chemical equilibrium holds when the velocity 
of the phase interface is much smaller than the diffusion velocity of baryon number 
and the length scales of the disconnected hadron or quark regions are much smaller 
than the baryon number diffusion length-scale during the transition. This is the 
case, for example, soon after bubble nucleation, when most of the bubbles present 
have radii around the critical value 

Rc=^^. (6) 

Ph-Pq 

where a is the surface tension associated with the phase interface. 

After nucleation, the bubble starts to expand (the low temperature phase is 
thermodynamically favoured) and its surface accelerates until it reaches a steady 
state velocity corresponding to self similar growth for a weak deflagration bubble 
||l6| . (We are working here within the scenario of small supercooling, in which weak 
deflagrations are the relevant mechanism for the motion of the phase interface, and 
will assume that the separation between the nucleation sites is not very small.) The 
expansion produces a deviation away from chemical equilibrium and the baryon 
number density can be increased in the vicinity of the interface piling up mostly 
on the quark side but also to some extent in the hadron phase as a result of the 
increased flux coming from the quark medium. The situation in the bulk of the 
two phases is not significantly modified, however. Diffusion tends to impede the 
accumulation of baryon number and, when the self similar growth stage is reached 
for the hydrodynamic variables, the baryon number density attains a stationary 
profile. (Note that the intrinsic length scale set by the diffusion coefficient prevents 
the baryon number density following a self similar solution.) Most of the excess 
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baryon number on the quark side of the interface is then within a layer of thickness 
with the value of the baryon number density there being joined to the background 
value in the bulk of the quark phase via a decaying profile. A rough estimate for 
is given by 

r.^-, (7) 

^/ 

where D is the baryon number diffusion coefficient and vj is the steady state velocity 
of the front as measured from the centre of the bubble. 

The self similar growth is first broken by the coupling between the radiation 
fluid and the standard fluids when the bubble has reached dimensions of the order 
of A [0. After the coupling however, the similarity solution for bubble growth is 



restored and it ultimately ceases when the spherical compression waves preceding 
the deflagration fronts of neighbouring bubbles start to interact. (This occurs well 
before the bubble surfaces themselves come into contact.) In principle, for a perfect 
fluid, the compression waves would be fronted by shocks but it should be noted that 
the predicted amplitude of such shocks is negligibly small for the case of spheri- 



cal bubbles and small supercooling which we are considering here [0, |T8|, 0, p!9| . 
However, independently of whether there is any significant shock or not, the kinetic 
energy of the ordered motion will be progressively converted into internal energy via 
compression, once adjacent bubbles have started to interact, and the temperatures 
of both phases will be raised near to the critical temperature Tc. At this stage, 
the bubbles grow much more slowly, on a time scale which is essentially set by the 
expansion of the Universe leading to cooling which allows for continuation of the 
transition. In this slow growth stage, global chemical equilibrium is restored and 
the baryon number which has accumulated near to the phase interface can be re- 
distributed into the bulk of the two phases. When the surfaces of adjacent bubbles 
meet, they coalesce to form larger bubbles thus minimizing the total surface energy. 
This coalescence gives rise to disconnected quark regions which then proceed to 
evaporate, tending to become spherical under the action of surface tension. 
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The hydrodynamical evolution of isolated evaporating quark drops has been 
investigated in detail in a number of studies [2^, 0, 0, |19| and can be essentially 
summarized as consisting of a self similar stage followed by one in which long range 
energy and momentum transfer takes place and then a final decay which may be 
dominated by surface tension effects and become increasingly rapid. 

During the first stage, quark drops evaporate converting quarks into hadrons at 
the rate necessary to keep the internal compression constant and uniform. In a very 
pictorial description, each quark drop can be viewed as behaving like a shrinking 
"leaky balloon" which is ejecting material at the rate necessary to avoid producing 
any compression. If baryon number were entirely carried along with the hydrody- 
namical flow and no suppression of baryon number flux occurred at the interface, 
then this stage of the evaporation would not produce any increases in the baryon 
number density. If, however, there is some flux suppression, this stage will still con- 
sist of a self similar evolution for all of the hydrodynamical variables apart from the 
baryon number which will accumulate ahead of the interface with diffusion to either 
side of it in a way similar to that discussed above for a growing hadron bubble. 

When the quark drop reaches dimensions comparable with the mean free path 
for the particles of the radiation fluid, the process of entropy extraction breaks 
the self similarity of the evaporation and produces a sharp increase in the fluid 
compression (of both the quark and hadron phases) []1T| and also in the baryon 
number density. It should be stressed that this increase in baryon number density is 
distinct from the one produced by flux suppression and is not necessarily localized 
at the phase interface but is rather extended over the whole quark region. This is 
because its origin is not hydrodjTiamical but is the consequence of the long range 
energy and momentum transfer via the radiation fluid particles. 

As pointed out in |14|, the very final stages of the drop evaporation (for 
Rg < 10^ fm) may be dominated by the surface tension a which could give rise 
to an increasingly rapid evaporation and a consequent increase of both the com- 
pression and the baryon number density. However, recent lattice gauge calcu- 
lations seem to indicate rather small values for the surface tension coefficient 
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(Jo = cr /T^ ^ 0.01 — 0.02 |jT8[ and, if this is correct, then surface tension would 
play only a minor role for the increase of baryon number density. 

In the next Section, we introduce the set of relativistic hydrodynamical equa- 
tions which we have used for describing the dynamics of single isolated evaporating 
quark drops and from whose numerical solution we have derived quantitative results 
for the profile and amplitude of the baryon number inhomogeneity produced at the 
end of the quark-hadron transition. 



III. Relativistic hydrodynamical equations 

We here briefly review the relativistic equations which we have used for calcu- 
lating baryon number concentration at the end of the quark-hadron transition and 



refer the reader to previous papers [12, 13, 14, 11| for further details. The equations 



can be divided into three main groups, related to: i) the dynamics of the spherical 
drop, ii) radiative transfer, in) the behaviour of the diffusing "baryon number fluid". 
We use a Lagrangian spherically symmetric coordinate system comoving with the 
fluid and having its origin at the centre of the drop. It is appropriate to write the 
metric line element as 



where /i is the comoving radial coordinate and R is the associated Eulerian coor- 
dinate. Local conservation of energy and momentum for the combined fluids (the 
standard fluid together with the radiation fluid), as measured in the local rest frames 
of the standard fluid, can be expressed as the vanishing of the projections of the 
combined stress-energy tensor parallel and perpendicular to the direction of the fluid 
four-velocity u. Using this, the continuity equation and the Einstein equations, the 
full set of equations for determining the behaviour of the standard fluids can be 
written as: 
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u t = —a 



e,t = wp,t - aso, (10) 



aw pw 



(12) 



M,^ = 4rfi?,^(^e + u7o + ^wi), (13) 
^ n=-i?,t, (14) 



47ri?V ' a 

r = -fi,= (i + .'-— ) . (15) 

where w = {e + p)/ p is the specific enthalpy of the standard fluids, u is the 
radial component of the fluid four velocity in the associated Eulerian frame, F is 
the general relativistic analogue of the Lorentz factor and M a generalized mass 
function. The quantity p is the compression factor which expresses the variation in 
proper volume of comoving elements of the standard fluids with respect to a fiducial 
value. (For a classical standard fluid, p can be taken to represent the rest mass 
density). The interaction between the standard fluids and the radiation fluid enters 
through the source functions Sq and Si and through the radiation contributions 
to the gravitational terms via the radiation energy density Wq, energy flux Wi and 
anisotropy W2- The latter quantities are calculated by considering local conservation 
of energy and momentum of the radiation fluid alone, and this provides the following 
set of equations |T^ 
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(wo),t + -^(wi),^ + -(^-^ + 



2a / a 



ova 



R 



R 



b t Rt 



W2 



asQ, 
(16) 



a /I 
6 U 



4a M ^9/^^t^^,t\ 



b \ a 



+ 



R 



W2 



asi, 
(17) 



W2 = /e^o 



where (18) is a closure relation written in terms of a variable Eddington factor 
In our calculations, the source functions Sq, Si are expressed as 



So 



:i+«2) 

A 




Wo 



"Si = — r 



A 



(19) 



with Qr and Tp being the number of degrees of freedom and the temperature of the 
radiation fluid. The parameter a2 is an adjustable coefficient (ranging between zero 
and one) which takes account of the contribution of non-conservative scatterings to 
the energy source function (for the present calculations we take a2 = 1). We treat 
the phase interface as a discontinuity surface and track it continuously through the 
finite difference grid using junction conditions which are imposed across the interface 
and a characteristic form of the above hydrodynamic equations which is used for the 
flow regions adjacent to it. A detailed discussion of these extra equations and of the 



way in which they are obtained has been given in previous papers [p!2| , |T3|, |TT|] and 
we will not repeat it here but, rather, just list the expressions used for the junction 
conditions. The energy and momentum junction conditions for the standard fluids 
are expressed as 



[{e + p)abf = 0, 



(20) 
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2 [abdtK f J ah fR 



As junction conditions for the energy and momentum of the radiation fluid after the 
decouphng from the standard fluids has begun to take place ITTi we use: 



± 



while for the metric coefficients we have: 



(22) 
(23) 



[R]^ = 



(24) 



[au + bfisT]^ = 



(25) 



2 ,-,2i± 



. 



(26) 



Note that the use of brackets in equations (|20|)-(p6D follows the conventions: [A]^ = 
y4+ — A~ , {A}^ = A'^ + A~ , and that the superscripts ^ indicate quantities imme- 
diately ahead of and behind the interface. The interface location is denoted by fig 
with fig = d^s/dt and f = {a? — b'^fi'^^Y^'^. 

Finally, for a phase interface which behaves as a weak deflagration front, it 
is necessary to supply one further equation giving the rate at which material is 
converted from one phase to the other (or, alternatively, to specify the interface 
velocity). As discussed in previous papers (see, for example [0]) we use a formula 
giving the rate of conversion of quarks into hadrons in terms of a black-body model 
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for the phase interface, suitably corrected with an adjustable accommodation coef- 
ficient «! (0 < ai < 1) to take account of deviations away from an ideal black-body 
law. We then have 



awfis 



T 



{9h + gr 



TT 

30 



(27) 



For the calculations presented here, we have taken ai = 1, but an extend discussion 
of computations with smaller values of ai = 1 can be found in 111. 



IV. Baryon number flux and diffusion 

As mentioned in the Introduction, simple considerations seem to suggest that 
baryon number would not be carried along entirely together with the hydrodynami- 
cal flow but, rather, that it would tend to be "trapped" inside the high temperature 
phase because of flux suppression at the phase interface. In order to study the evo- 
lution of the baryon number distribution, we need to extend the multi-component 
fluid description discussed in Section III by adding to the treatment of the standard 
fluids and the radiation fluid, an additional "baryon number fluid". In particular, 
we need to introduce equations to describe the hydrodynamics of a fluid which has 
a suppressed flow at the phase interface and which diffuses relative to the standard 
fluids. 

The effects of this diffusion can be included by introducing a diffusive flux 
four-vector q into the baryon number continuity equation which becomes 

(n6M" + g")^„ = 0, (28) 

where is the baryon number density, u is still the four velocity of the standard 
fluid and Greek indices are taken to run from to 3. In the frame comoving with 
the standard fluid, = (1/a, 0, 0, 0) and q has only a spatial component [i.e. 
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UaQ°' = and = (0, q^, 0, 0)]. Using the property (V^"); a = iV—gV^), a/v^— ^ 
(where g is the determinant of the metric tensor) and the metric relations @, (|I5]) 
we can rewrite (EHI) as 



bR' 



dt 



\R^J OR 



aR^ djbT) djPnb) 
~r M OR 



R^ djPnb) 
br OR 



, (29) 



(see the Appendix for details) where the radial component of the diffusive flux is 
written as 



1 d{Dn,) 



1 d{Dnb) 



^ " {bvy dfi bT OR ■ ^ ' 

In the Newtonian limit, T = a = 1 and equation ( pOj) reduces to the standard 
Lagrangian diffusion equation in spherical coordinates, with D being the diffusion 
coefficient (note that here, d/dt is a Lagrangian time derivative for a location co- 
moving with the standard fluid). An attractive feature of Lagrangian schemes is 
that advection is treated exactly. Exploiting this, it is convenient to use for the 
finite-difference representation of equation (|29|), the simple and compact expression: 



47r(nfe^"+^ 



i+i/2 



/p2\n+l 
^ J+1/2 



i+i/2 



An 

pTl+l/2 

^ i+i/2 



i+i/2 



^ i+1/2 



'J+l/2 



(31) 



where the superscripts refer to the time level at which the quantity is calculated and 
the subscripts to the position in the spatial grid, AR^_^_^^2 = ^j+i ~ -^i 



($ 



D 



-D 



i+i/2 



{rib) 



i-i/2 



pn _ rpn 



(32) 



14 



is the diffusive flux of baryon number. A rougli estimate for the value of the diffusion 
coefficient D (which we take to be constant in time, uniform in space and the same 
in both phases) can be deduced with rather simple arguments if we rule out the 
possibility of the diffusion being turbulent (which seems likely to be correct 0). 
In this case, baryon number diffusion can be described as a simple Brownian motion 
of baryon number carriers having a mean free path of the order \free ^ |^ 
giving a microscopic diffusion coefficient D ~ 10^^ — 10 fm. 

Clearly, the baryon number within a given grid zone changes in time only if 
a baryon flux crosses its boundaries. While for the bulk of the two phases, the 
flux leading to the variation is just the diffusive flux (|32|) , a special treatment is 
necessary for the two parts of the grid zone containing the phase interface. At the 
interface, the diffusive flux is accompanied by a much larger flux which is related to 
the hydrodynamical flow <l>p of elements of the quark gluon plasma as it is converted 
to the hadron phase. As mentioned earlier, this flux could be subject to suppression 
processes at the interface but no exact expression is yet available for this. We 
therefore proceed here by defining a phenomenological expression for the net baryon 
number flux across the interface $b in terms of the hydrodynamical flux $p and 
a suitable "filter factor" F which expresses the ratio between the baryon number 
passing across the phase interface and the total baryon number incident on it. 

In principle, F could be expressed in terms of the probability of finding (from 
all of the quarks and antiquarks present) three quarks of the right types within 
a volume of 1 fm^ and in a time of 10~^^ s, in terms of the "transparency" of the 
phase interface to the passage of a baryon number carrier from the quark phase to the 
hadron phase T^q^h and of the corresponding probability S/^^g that a baryon hitting 
the phase boundary from the hadron phase is absorbed Moreover, referring to a 
situation in which chemical equilibrium holds, it would be possible to express T^g^h 
in terms of S/^^g and this would restrict the uncertainty to this latter quantity 
only. (Although convenient, this approach requires that the baryon transmission 
probability does not vary significantly when the equilibrium is broken.) 
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Unfortunately, no reliable value for the baryon transmission probability S/^^g is 
known at present and, worse than this, different approaches to the study of the rates 
of elementary processes taking place at the phase interface seem to result in quite 
different estimates of it (see P, |2l| and [^] for further references). In view of this 
uncertainty, we will treat the filter factor as essentially a free parameter, adopting 
the reference value F ~ 10^^ as estimated from the expressions presented by Fuller 
et al. for = 150 MeV and S/i^^ ~ 10^'^. (The value F = 1 corresponds to 
the case where the baryon number flow crosses the phase interface unimpeded and 
clearly represents an upper limit.) 

The flux of elements of the standard fluid across the interface can be evaluated 
by projecting the flux four-vector along the unit spacelike four-vector n normal to 
the timelike hypersurface describing the time evolution of the interface. We obtain 

where /is < 0, / = (a^ - b'^filY^'^ and 

= -(1, 0, 0, 0), n„ = ^(-As, 1, 0, 0) . (34) 

a J 

Neglecting diffusive contributions to the baryon number flux across the interface, 
we obtain 



and we assume continuity of this flux across the interface (i.e. = 0). 

In the next Section we will present results from numerical computations made 
using the set of equations which we have introduced here and illustrate the roles 
played in the segregation of baryon number by suppression mechanisms and by 
radiative transfer. 
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V. Numerical strategy and results 



The numerical approach implemented for the solution of the present set of 
equations (D-(|15D, (|16D-(|19D, (|0D-(||) and (H) is based on that used for our 
earlier computation of the evaporation of a quark drop including the effects of long 



range energy and momentum transfer [|TT] with an extension to include the parabolic 
diffusion equation (p9|). 

As is frequently the case for numerical computations in which a diffusion equa- 
tion needs to be solved, we are here faced with the problem of performing calcu- 
lations with a mixed set of equations and avoiding excessive restrictions placed on 
the time step by the von Neumann criterion for stability in solving the parabolic 
equation. Obeying this criterion, which is more stringent than the usual Courant 
one since it depends quadratically on the minimum grid spacing and is inversely 



proportional to the diffusion coefficient allows one to perform an explicit inte- 



gration of equation (|3l|) without any further constraint. However, for the present 
case, implementing the von Neumann condition produced computational times for 
each simulation which were not affordable. The reason for this is connected with 
the particular organization of our grid which has an exponentially increasing spacing 
in order to facilitate following the dynamics of the drop through several orders of 
magnitude change in the radius. To avoid this problem, we have implemented a 
standard "flux limiter" scheme in which a control is set on the diffusive flux in the 
quark phase preventing it from evacuating the baryon number content of any grid 
zone within a single time-step [ES]. The results obtained in this way were found to 
be in excellent agreement with ones from a comparison calculation without the flux 
limiter and using the von Neumann condition (but which required a computational 
time longer by a factor of twenty). 

Another concern in the present calculations has been that of preserving as 
closely as possible the overall conservation of baryon number. Clearly results for the 
final distribution of baryon number at the end of the transition will be worthless if the 



integration of the diffusion equation (31) is not accurate enough and is significantly 
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producing or destroying baryon number. This equation is nearly conservative in form 
but special attention needed to be paid to the calculation of interpolated values for 
a and F (for which we used function fitting procedures) and to the implementation 
of the regridding procedure pi| , [TT| . Having done this, our computations preserve 
the total baryon number in the grid to an accuracy of a few parts in 10^ for runs of 
about 10^ time steps. 

In the following, we present results from computations which have been per- 
formed evolving from initial data given by the self similar solutions derived in |T^ 
for a spherical quark drop with initial dimensions Rsfi = 10^ fm, at an initial tem- 
perature Tq = Tq/Tc = 0.998, surrounded by a hadron plasma at temperature 
Th = Th/Tc = 0.990. The similarity solutions provide suitable initial conditions for 
all of the variables apart from the baryon number which does not necessarily follow 
the bulk hydrodynamical flow. We decided to start with the baryon number density 
in each phase being uniform and given by the expressions (0) and (^) presented in 
Section II which correspond to conditions of chemical equilibrium. This is very ap- 
proximate because global chemical equilibrium does not apply for a situation with a 
moving interface such as the one which we are considering but, nevertheless, initial 
conditions imposed in this way are sufficiently good to allow the solution to relax 
rapidly to a consistent one. 

We have examined the effects on the overall solution of varying the values of 
important input parameters and we will be discussing this in detail but first we 
concentrate on the results obtained for a set of fiducial parameter values. This 
"standard" run follows the evaporation of a quark drop with surface tension pa- 
rameter (To = 0.01, filter factor F = 0.3 and diffusion coefficient D = 1 fm. The 
decoupling radius Rd (i.e. the drop radius at which the decoupling between the 
standard fluids and the radiation fluid is allowed to take place |Tl[]) is set to be 10"^ 
fm, which corresponds to the average mean free path of the electromagnetically- 
interacting particles. We do not consider here a decoupling with neutrinos which, 
however, would follow a similar hydrodynamical behaviour except for the different 
number of degrees of freedom involved [jlTI . 
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Figure 1. Time evolution of the most significant hydrodynamical vari- 
ables. Each curve refers to a given time in the solution, with the quark 
phase being to the left of the vertical discontinuity and the initial con- 
ditions being indicated with the dashed curves. Starting from the upper 
left-hand window and proceeding clockwise, the frames show: the radial 
component of the fluid four-velocity in the Eulerian frame u, the energy 
densities of the standard fluids e and of the radiation fluid u'o and the 
compression factor p. The decoupling between the radiation fluid and 
the standard fluids is allowed to start at Rg = fm 
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Figure 2. Time evolution of the logarithm of the baryon number density. 
The continuous curves refer to successive stages in the evolution, with 
the initial conditions being indicated by the dashed curve. 
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Figure 3. Three-dimensional plot of the profile of the log of the baryon 
number density at the end of the computation, when the quark drop 
had a radius Rg = 2 fm. 
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This decoupling would take place at a scale X,y ~ 10^^ fm but it is not currently clear 
whether spherically isolated quark drops would have been present at that scale. 

Figure 1 shows the time evolution of the most important hydrodynamical vari- 
ables: the radial component of the fluid four-velocity u, the energy density of the 
standard fluids e, the compression factor p and the energy density of the radiation 
fluid Wq. Although similar curves have already been presented and discussed in 
TT| , we consider it useful to show them here as they provide the hydrodynamical 
background for our discussion of baryon number segregation. 

The main features of the hydrodynamical solution can be summarized as fol- 
lows: i) the presence of a self similar solution until the radiation decoupling starts 
at Rs = Rd] ii) the recovering of an "almost" self similar solution for the energy 
density e and the fluid velocity u after decoupling; in) the significant increase of the 
compression factor in both phases caused by the extraction of entropy by the radi- 
ation fluid; iv) the smoothing out of the step in the profile of the radiation energy 
density when decoupling occurs. 

In Figure 2 we show the time evolution of the baryon number density which, 
for convenience, we can divide into three main stages. In the first of these, baryon 
number density deviates from the imposed initial conditions and settles into a so- 
lution which is consistent with the motion of the interface. As mentioned above, 
the initial conditions adopted assume chemical equilibrium to hold across the phase 
interface. However, when the interface is in motion (as it is for our initial condi- 
tions), equilibrium cannot hold and the baryon number density in the hadron phase 
has to adjust itself in order to reach a self consistent solution. By the time the drop 
radius has decreased from 10^ fm to approximately 2.5 x 10^ fm, the new solution 
has reached its steady state configuration. 

The next stage of the evolution is characterized by a "snowplow" behaviour 
of the phase interface which accumulates baryon number as it moves into the quark 
medium. This piling-up of baryon number on the quark side is counteracted by the 
effects of diffusion which tends to smear out the build-up, producing a rapidly decay- 
ing profile (a similar behaviour was also suggested by Kurki-Suonio [^). During this 
stage, which corresponds to a self similar solution for the hydrodynamical variables 
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presented in Fig. 1, the counteracting mechanisms of baryon number flux suppres- 
sion and of baryon diffusion, suitably regulate the rate of baryon number passage 
across the phase interface until a stationary solution is reached for the baryon num- 
ber density. Note that the length scale is too small to be resolved on the grid at 
this stage and this explains why, in Fig. 2, the decaying profile appears to have a 
constant logarithmic width (connected with the grid structure). 

The self similar solution for the hydrodynamical variables can be maintained 
only as long as the quark drop evaporation is effectively scale free and no other 
process is intervening. Equivalently, the stationary solution for baryon number can 
be maintained only as long as the underlying hydrodynamical behaviour remains 
self similar and the tail of the diffusing profile does not interact with the centre of 
the drop (the time when this happens is dependent on the magnitude of D) . In the 
present situation there are two possible length scales that when reached by the drop 
radius would produce a deviation away from the hydrodynamic similarity solution. 
These are the mean free path of the radiation fluid particles (at which scale the 
radiative transfer is most effective) and the length scale set by the surface tension 
(which is much smaller). 

For our fiducial parameter values, the radiative transfer occurs before the drop 
can know about its centre via diffusion and, as a consequence, the following stage 
of the evolution in Fig. 2 starts at Rg — 10^ fm. When we compare the effects of 
the entropy extraction via radiative transfer on the compression factor (lower-right 
window in Fig. 1) and on the baryon number density, we can find analogies and 
differences. The first common feature is the increase of p (and of rib) o-lso in the 
low temperature phase for Rs < 10^ fm. However, while the radius at which the 
compression increase in the hadron phase begins does not change in time, this is 
not the case for the radius at which the increase in baryon number density begins. 
The latter slowly shifts towards larger values under the action of baryon diffusion. 



22 



Figure 3 shows a three-dimensional representation of the profile of baryon number 
density at the end of the computation when the quark drop has radius Rg = 2 fm. i 

Another difference between the evolution of the compression factor and that 
of the baryon number density is the larger increase in the latter (about 9 orders 
of magnitude) produced after radiation decoupling. This is due to the fact that 
baryon number density can be increased not only by the entropy extraction via 
radiative transfer but also by flux suppression at the interface and the latter is a 
much more effective mechanism. Note that the profile of baryon number density 
shown in Fig. 3 corresponds to material which has non-uniform specific entropy 
(as a result of radiative transfer but which is very nearly in thermal and 

mechanical equilibrium in each phase |P, |10| . 

When looking at the profiles in Figures 2 and 3 it is not possible to distinguish 
the different contributions coming from the radiative transfer and from the flux 
suppression. However, using the mean free path of the radiation particles A as a 
free parameter, it is possible to regulate the action of the radiative transfer. Taking 
A = corresponds to a situation in which the decoupling never takes place at all. 
Figure 4 shows the final profile of the baryon number density when the quark drop 
has radius i?^ = 2 fm as calculated for progressively smaller values of the mean free 
path, which is equal to A = 10^ fm for the heavy continuous curve, A = 10^ fm for 
the dashed curve and to zero for the dotted one. 

The effects of radiative transfer can be clearly seen by comparing the dotted 
curve with the heavy continuous curve and amount to a relative increase in the 
baryon number density of around one order of magnitude in the inner part of the 
profile. Note that the use of a logarithmic scale can be misleading and seem to 
suggest that total baryon number in the different plots is not the same. However, 
when comparing profiles for different values of A, one should pay particular attention 
to differences in the outer regions which are not very noticeable but which make a 

■^It is important to stress that, we here stop our computations at Rg = 2 fm because 
already at this stage the whole hydrodynamical approach described in Sections I-III is 
ceasing to be valid. It is no longer sensible to use a fluid description for a quark drop 
containing only very few quarks. 
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considerable contribution. In order to emphasize these differences, a magnification 
of the final profiles is shown in the small diagram at the top right of Figure 4. 
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Figure 4. Final profile of the baryon number density when the quark 
drop has radius i?s = 2 fm as calculated for different values of the mean 
free path A. The curves correspond to A = 10*^ fm (heavy continuous 
curve), 10^ fm (dashed curve) and (dotted curve). The small diagram 
at the top right shows a magnification of the differences between the 
curves. 

Figure 5 shows results from calculations performed with different values of the 
filter factor, showing the final profiles of baryon number density when the drop has 
radius 2 fm. Different curves refer to values of the filter factor F ranging from 0.2 
to 1, with the heavy continuous curve corresponding to the preferred value F — 0.3. 

The behaviour of the solution under variation of the filter factor F is straight- 
forward to understand: a more "transparent" phase interface (i.e. with a higher 
value of F) leads to less accumulation of baryon number inside the drop and pro- 
duces a smaller final increase in the baryon number density. However, because of 
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the nonlinearity of the process, a less transparent interface actually gives a larger 
net baryon number flux as a consequence of the increased baryon number density in 
the quark phase [see equation (|35D]. This explains why the baryon number density 
is larger also outside the quark drop. Note that also in this figure, the logarith- 
mic scaling could be confusing and suggest that the total baryon number is not the 
same for the different curves. However, once again, as shown in the small magnified 
diagram, the large differences which appear for Rs < 10'^ fm are compensated by 
changes further out which are less easy to see. 
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Figure 5. Final profile of the baryon number density when the quark 
drop has radius Rg = 2 fm as calculated for different values of the filter 
factor at the interface F. The curves correspond to F = 0.2 (dashed 
curve), F = 0.3 (heavy continuous curve), F = 0.4, 0.5, 0.6, 0.7 (stan- 
dard continuous curves) and F = 1 (dotted curve). The small diagram 
at the top right shows a magnification of the differences between the 
most important curves. 

A number of comments should be made about the final baryon number density 
reached inside the quark phase for different values of the filter factor. As shown by 
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the curve corresponding to F = 0.2, the baryon number density in the quark phase 
could well reach values above that for nuclear matter and this could suggest that a 
transition to strange quark nuggets might possibly occur at some stage during the 
quark drop evaporation. This is a very speculative idea and we limit ourselves to 
pointing out that, as far as baryon number segregation is concerned, conditions for 
the creation of quark nuggets are not difficult to reach if a strong suppression of 
baryon number flux takes place at the interface. Another important point to notice 
is that with the present choice of parameter values, only a very small fraction of the 
initial total baryon number will remain in the high density regions after drops have 
evaporated away. A final comment should be made about the properties of the final 
baryon number density profile obtained for F = 1, which corresponds to a situation 
in which baryon number carriers are freely streaming across the interface and is 
the result of radiation decouphng only. Although this case may well not refer to a 
realistic scenario, it is instructive as it shows the role of the baryon flux suppression 
in the production of baryon number density peaks. It is interesting to compare this 
with the A = curve of Fig. 4 which corresponds to radiative transfer not being 
operative. 

We conclude our analysis of the numerical results by presenting the final baryon 
number density profiles obtained after varying the diffusion coefficient while main- 
taining the other parameters with their standard values. Figure 6 shows curves 
calculated for values of D ranging between zero (dotted curve) and 10 and with the 
heavy continuous curve corresponding to the preferred value D — 1. 

The curves demonstrate the action of diffusion which prevents the accumula- 
tion of baryon number in a narrow shell ahead of the phase interface (as seen for 
D — 0) and spreads out the accumulation of baryon number produced in the hadron 
phase near to the drop surface. Larger values of the diffusion coefficient give rise to 
larger final values of the baryon number density in the quark phase, but the results 
show that the high number density regions are unlikely to contain a large fraction of 
the total baryon number if only microscopic random diffusion is taken into account 
(which leads to relatively small values for D) . 
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Figure 6. Final profile of the baryon number density when the quark 
drop has radius Rg = 2 as calculated for different values of the 
diffusion coefficient D. The curves correspond to -D = 10, 1 (the heavy 
continuous line), 0.1 and (the dotted curve). 

Stronger diffusion is more efficient in preventing accumulation of baryon num- 
ber adjacent to the quark side of the interface and this, in turn, tends to reduce 
the outward flux of baryon number into the hadron phase. The effective role played 
by diffusion is then that of trapping baryon number within the high density region 
and so it is no surprise that higher final values of nj, should result for larger values 
oi D. It is worth underlining that the values adopted for our simulations refer to a 
microscopic modeling of diffusive processes in which turbulent motion is neglected. 
It is possible that non radial fiuid motions could produce a larger effective diffusion 
coefficient, but not enough is yet known about this. 
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VI. Conclusions 



We have presented here results from numerical computations of the final stages 
of a first order cosmological quark-hadron phase transition, aiming to estimate the 
baryon number density profile produced. For doing this, we have used mathemat- 
ical and numerical techniques previously developed for following the evaporation 
of spherical quark drops together with a separate treatment of a distinct fluid of 
baryon number carriers. The effects of baryon number flux suppression at the phase 
interface and of baryon number diffusion have been included within the overall rel- 
ativistic hydrodynamical treatment. The results obtained have shown that baryon 
number segregation is an inevitable result of a first order transition and illustrate 
the different roles played in baryon number segregation by flux suppression and by 
radiative transfer. Final baryon number density profiles have been computed which 
exhibit relative increases of several orders of magnitude in both phases, with an 
overdense spherical region of radius R ^ X ^ 10^ fm being left behind in thermal 
equilibrium in the hadron phase. This region, which probably represents the most 
important relic of the transition, will be subject to further subsequent diffusion. 

As pointed out by Ignatius et al. |]18|, a number of conditions would need to be 
met if baryon number segregation is to have any significant effect on nucleosynthe- 
sis. One of these conditions (that the final baryon number density contrast should 
be larger than ~ 10) is certainly satisfied by the results from our present study. 
However, our results also indicate that the condition for the high density regions 
to contain most of the total baryon number is extremely difficult to satisfy unless 
turbulent diffusion might be operative. A third necessary condition requires that the 
mean separation scale In between two baryon number density peaks of the type com- 
puted here should be rather large. It would be necessary to have /„ > 1 m ^ 10^"^ t^, 
where tn is the horizon scale at the time of the transition. As mentioned in the In- 
troduction, this is a fundamental quantity for the whole description of the phase 
transition, without which any conclusion about the role played by baryon number 
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inhomogeneities on the subsequent nucleosynthesis remains extremely unclear. Con- 
siderable efforts have been devoted to determining it (see Christiansen and Madsen 
for a recent discussion) but the value of this scale is still uncertain and further 



developments need to be awaited. 

The hydrodynamical scenario for baryon number segregation discussed here, 
could represent a promising starting point for studying the possible production of 
seed magnetic fields during the quark-hadron phase transition. Cheng and Olinto 
P3| pointed out that because of the different fractional electric charge carried by up 
and down quarks, charge separation could be produced together with the concentra- 
tion of baryon number in the high temperature phase. In their paper, they focussed 
attention on the slow growth stage of the phase transition during which chemical 
equilibrium holds and the baryon number contrast can be estimated analytically. 
However, this same mechanism might be at work also during the final stages of 
the transition, for which the present calculations provide more precise quantitative 
estimates of the dynamical motions of charges and of their density. 

Because of the much more extreme baryon number concentrations arising in 
the present context, the magnetic field could be several orders of magnitude larger 



than that discussed in ||2^ and might possibly reach the equipartition value. This 
requires further investigation and will be a focus of future work. 
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Appendix 



In this Appendix, we briefly sketch some intermediate steps leading to equation 
(p9|), the relativistic diffusion equation for baryon number. Using the standard 
property of four-divergences mentioned in the text, equation (E8|) can be written as 



bR' 



{nb),t + nb — + 



-ah 



\ a J 



b R 

The T° component of the Einstein field equation provides the useful relation 
bt R 



(36) 



^■, fit 



R 



aR 



a 



(37) 



where G <^ 1 is the gravitational constant and we drop the term containing it for 
the present purposes. Making use of the identity 



b, fi _ R, fxfj. ^ 
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and of the definition (|30|) , it is then possible to rewrite equation ( p6l) in the final 
form (|29|). 
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